map-examples
Steve Simon
December 29, 2016
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("ggmap")## Warning: package 'ggmap' was built under R version 3.3.2## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.3.2library("lubridate")## Warning: package 'lubridate' was built under R version 3.3.2##
## Attaching package: 'lubridate'## The following object is masked from 'package:base':
##
## datelibrary("magrittr")##
## Attaching package: 'magrittr'## The following object is masked from 'package:ggmap':
##
## insetlibrary("tidyr")## Warning: package 'tidyr' was built under R version 3.3.2##
## Attaching package: 'tidyr'## The following object is masked from 'package:magrittr':
##
## extractlibrary("tibble")
r <- 6371
p <- 2*pi*r / 360
a <- cos(39*pi/180)Latitude starts at 0 at the equator, goes up to 180 degrees at the North Pole, and down to -180 degrees at the South Pole. Longitude stats at 0 at the Greenwich Meridian, goes up to 180 as you head west and down to -180 as you head east. The radius of the earth is roughly 6371. Multiply this by 2 pi and divide by 360 to get 111.2 kilometers per degree of latitude. Longitude is a bit trickier. At the equator, a degree of longitude is the same as a degree of latitude. But things shrink as you move towards the poles. You have to adjust for this, and the adjustment factor equals the cosine of the latitude. So, for example, Kansas City is roughly at 39 degrees of latitude, so the adjustment factor is 0.777, making a degree of longitude only 86.4 kilometers.
This is a reasonable approximation, but you can adjust for the fact that the earth is not a perfect sphere (it is slightly flatter at the poles than at the equator). You can also adjust for the curvature of the earth if the distances are far enough apart. The R bloggers site has a [https://www.r-bloggers.com/great-circle-distance-calculations-in-r/ nice summary] of the more accurate formulas.
One last thing: I usually get this wrong, but when you draw a plot, longitude goes on the x-axis and latitude goes on the y-axis.
Example 1. A gpx file.
I run a couple of miles every second or third day and participate in a few timed races (mostly 5K and 4 mile events) from time to time. I record my runs using an iPhone app, MotionX-GPS. It produces an xml file that includes geographic positions and time throughout the run. Although MotionX-GPS produces nice plots of my runs, I wanted to produce something a bit more detailed using R. Here’s an example using the data from a 5K race on January 1, 2017.
The file Track 502.gpx should be available at my github site.
The “interesting” lines in this file look something like
“
or
“
or
“”
There are other lines in the file which you can ignore safely. The first step is to pull out the relevant pieces of information amid all the xml code. I use regular expressions to do this.
"Track 502.gpx" %>% read.delim(header=FALSE, as.is=TRUE, sep="~") -> gpx_lines
gpx <- NULL
gpx_lines$V1 %>%
grep("<time>", ., value=TRUE) %>%
sub(".+T", "", .) %>%
sub("Z</time>", "", .) %>%
as.numeric %>%
as_tibble %>%
bind_cols(gpx) %>%
rename(timd=value) -> gpx## Warning in function_list[[i]](value): NAs introduced by coerciongpx_lines$V1 %>%
grep("trkpt lat=", ., value=TRUE) %>%
sub("<trkpt lat=", "", .) %>%
sub(" lon=.+", "", .) %>%
as.numeric %>%
as_tibble %>%
bind_cols(gpx) %>%
rename(lat=value) -> gpx
gpx_lines$V1 %>%
grep("trkpt lat=", ., value=TRUE) %>%
sub(".+ lon=", "", .) %>%
sub(">", "", .) %>%
as.numeric %>%
as_tibble %>%
bind_cols(gpx) %>%
rename(lon=value) -> gpx
head(gpx)## # A tibble: 6 × 3
## lon lat timd
## <dbl> <dbl> <dbl>
## 1 -94.66584 38.87889 NA
## 2 -94.66583 38.87888 NA
## 3 -94.66579 38.87878 NA
## 4 -94.66574 38.87869 NA
## 5 -94.66568 38.87859 NA
## 6 -94.66562 38.87850 NAThe ggmap library has a “one stop shopping” function, qmplot (similar to qplot in ggplot2) that will automate the process of displaying geographic co-ordinates on a graph.
qmplot(lon, lat, data=gpx)## Using zoom = 17...## Map from URL : http://tile.stamen.com/toner-lite/17/31068/50147.png## Map from URL : http://tile.stamen.com/toner-lite/17/31069/50147.png## Map from URL : http://tile.stamen.com/toner-lite/17/31070/50147.png## Map from URL : http://tile.stamen.com/toner-lite/17/31071/50147.png## Map from URL : http://tile.stamen.com/toner-lite/17/31072/50147.png## Map from URL : http://tile.stamen.com/toner-lite/17/31068/50148.png## Map from URL : http://tile.stamen.com/toner-lite/17/31069/50148.png## Map from URL : http://tile.stamen.com/toner-lite/17/31070/50148.png## Map from URL : http://tile.stamen.com/toner-lite/17/31071/50148.png## Map from URL : http://tile.stamen.com/toner-lite/17/31072/50148.png## Map from URL : http://tile.stamen.com/toner-lite/17/31068/50149.png## Map from URL : http://tile.stamen.com/toner-lite/17/31069/50149.png## Map from URL : http://tile.stamen.com/toner-lite/17/31070/50149.png## Map from URL : http://tile.stamen.com/toner-lite/17/31071/50149.png## Map from URL : http://tile.stamen.com/toner-lite/17/31072/50149.png## Map from URL : http://tile.stamen.com/toner-lite/17/31068/50150.png## Map from URL : http://tile.stamen.com/toner-lite/17/31069/50150.png## Map from URL : http://tile.stamen.com/toner-lite/17/31070/50150.png## Map from URL : http://tile.stamen.com/toner-lite/17/31071/50150.png## Map from URL : http://tile.stamen.com/toner-lite/17/31072/50150.png## Map from URL : http://tile.stamen.com/toner-lite/17/31068/50151.png## Map from URL : http://tile.stamen.com/toner-lite/17/31069/50151.png## Map from URL : http://tile.stamen.com/toner-lite/17/31070/50151.png## Map from URL : http://tile.stamen.com/toner-lite/17/31071/50151.png## Map from URL : http://tile.stamen.com/toner-lite/17/31072/50151.png## Map from URL : http://tile.stamen.com/toner-lite/17/31068/50152.png## Map from URL : http://tile.stamen.com/toner-lite/17/31069/50152.png## Map from URL : http://tile.stamen.com/toner-lite/17/31070/50152.png## Map from URL : http://tile.stamen.com/toner-lite/17/31071/50152.png## Map from URL : http://tile.stamen.com/toner-lite/17/31072/50152.png## Map from URL : http://tile.stamen.com/toner-lite/17/31068/50153.png## Map from URL : http://tile.stamen.com/toner-lite/17/31069/50153.png## Map from URL : http://tile.stamen.com/toner-lite/17/31070/50153.png## Map from URL : http://tile.stamen.com/toner-lite/17/31071/50153.png## Map from URL : http://tile.stamen.com/toner-lite/17/31072/50153.png## Map from URL : http://tile.stamen.com/toner-lite/17/31068/50154.png## Map from URL : http://tile.stamen.com/toner-lite/17/31069/50154.png## Map from URL : http://tile.stamen.com/toner-lite/17/31070/50154.png## Map from URL : http://tile.stamen.com/toner-lite/17/31071/50154.png## Map from URL : http://tile.stamen.com/toner-lite/17/31072/50154.png## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead
You will see some additional examples using this data a bit later.
Example 2. KC 311 calls
The city government of Kansas City, Missouri has an open data initiative and one of their more interesting data sets is information from their 311 line. The 311 line is a phone number you can dial locally that will allow you to contact the city government to register a complaint or concern.
Please download the 311 center center service requests file from
https://data.kcmo.org/browse?category=311
fn <- "311_Call_Center_Service_Requests.csv"
cc <- read.csv(file=fn, header=TRUE, as.is=TRUE)
names(cc) %<>% tolower
str(cc)This file does not need any serious data manipulation. Here’s a quick default plot of the phyiscal location of 100 of the phone calls.
cc_sub <- cc[1:100,]
qmplot(longitude, latitude, data=cc_sub)Example 1, continued.
You can get and store a map that covers the longitude and latitudue values in the GPX file. This allows you a bit more flexibility in how you plot data on the map.
The bounding box is a range that helps insure that every longitude and latitude value falls inside the map.
bb <- make_bbox(lon, lat, data=gpx)
print(bb)## left bottom right top
## -94.66828 38.86852 -94.65759 38.88440The zoom value is a number between 1 and 22. A large number (like 3) would show an entire continent. A medium number (like 10) would show an entire city. A small number (like 17) would show a few city blocks.
"2401 Gillham Road Kansas City MO 64108" %>%
geocode %>%
get_map(zoom=17) %>%
ggmap## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=2401%20Gillham%20Road%20Kansas%20City%20MO%2064108&sensor=false## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.083761,-94.576988&zoom=17&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
"2401 Gillham Road Kansas City MO 64108" %>%
geocode %>%
get_map(zoom=10) %>%
ggmap## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=2401%20Gillham%20Road%20Kansas%20City%20MO%2064108&sensor=false## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.083761,-94.576988&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
"2401 Gillham Road Kansas City MO 64108" %>%
geocode %>%
get_map(zoom=4) %>%
ggmap## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=2401%20Gillham%20Road%20Kansas%20City%20MO%2064108&sensor=false## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.083761,-94.576988&zoom=4&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
You can ask for a value of zoom that is just the right size for your bounding box, though, for reasons that elude me, I have needed to adjust the zoom level downward by two levels.
zm <- calc_zoom(bb, adjust=-2L)
print(zm)## [1] 15Once you have a bounding box and a zoom level, you can get a map that will cover all of your data values.
Note the stucture of the objects created. The get_map function creates, by default, a 1280 by 1280 matrix of character strings that represent rgb color values in hexadecimal. Note the attributes that are stored with this matrix.
The ggmap function creates an object that is mostly a bunch of functions for displaying things. This fits within the framework of the ggplot2 library. You can display the map by itself just by typing the name of the object, or you can elements to the map with functions like geom_point or geom_text.
mp <- get_map(bb, zoom=zm)## Warning: bounding box given to google - spatial extent only approximate.## converting bounding box to center/zoom specification. (experimental)## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=38.876458,-94.662933&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=falsestr(mp)## chr [1:1280, 1:1280] "#FEFEFE" "#FEFEFE" "#FEFEFE" "#FEFEFE" ...
## - attr(*, "class")= chr [1:2] "ggmap" "raster"
## - attr(*, "bb")='data.frame': 1 obs. of 4 variables:
## ..$ ll.lat: num 38.9
## ..$ ll.lon: num -94.7
## ..$ ur.lat: num 38.9
## ..$ ur.lon: num -94.6
## - attr(*, "source")= chr "google"
## - attr(*, "maptype")= chr "terrain"
## - attr(*, "zoom")= num 15im <- ggmap(mp)
str(im)## List of 9
## $ data :'data.frame': 4 obs. of 2 variables:
## ..$ lon: num [1:4] -94.7 -94.6 -94.7 -94.6
## ..$ lat: num [1:4] 38.9 38.9 38.9 38.9
## ..- attr(*, "out.attrs")=List of 2
## .. ..$ dim : Named int [1:2] 2 2
## .. .. ..- attr(*, "names")= chr [1:2] "lon" "lat"
## .. ..$ dimnames:List of 2
## .. .. ..$ lon: chr [1:2] "lon=-94.67664" "lon=-94.64918"
## .. .. ..$ lat: chr [1:2] "lat=38.86575" "lat=38.88713"
## $ layers :List of 3
## ..$ :Classes 'LayerInstance', 'Layer', 'ggproto' <ggproto object: Class LayerInstance, Layer>
## aes_params: list
## compute_aesthetics: function
## compute_geom_1: function
## compute_geom_2: function
## compute_position: function
## compute_statistic: function
## data: waiver
## draw_geom: function
## finish_statistics: function
## geom: <ggproto object: Class GeomBlank, Geom>
## aesthetics: function
## default_aes: uneval
## draw_group: function
## draw_key: function
## draw_layer: function
## draw_panel: function
## extra_params: na.rm
## handle_na: function
## non_missing_aes:
## optional_aes:
## parameters: function
## required_aes:
## setup_data: function
## use_defaults: function
## super: <ggproto object: Class Geom>
## geom_params: list
## inherit.aes: TRUE
## layer_data: function
## map_statistic: function
## mapping: NULL
## position: <ggproto object: Class PositionIdentity, Position>
## compute_layer: function
## compute_panel: function
## required_aes:
## setup_data: function
## setup_params: function
## super: <ggproto object: Class Position>
## print: function
## show.legend: NA
## stat: <ggproto object: Class StatIdentity, Stat>
## aesthetics: function
## compute_group: function
## compute_layer: function
## compute_panel: function
## default_aes: uneval
## extra_params: na.rm
## finish_layer: function
## non_missing_aes:
## parameters: function
## required_aes:
## retransform: TRUE
## setup_data: function
## setup_params: function
## super: <ggproto object: Class Stat>
## stat_params: list
## subset: NULL
## super: <ggproto object: Class Layer>
## ..$ :Classes 'LayerInstance', 'Layer', 'ggproto' <ggproto object: Class LayerInstance, Layer>
## aes_params: list
## compute_aesthetics: function
## compute_geom_1: function
## compute_geom_2: function
## compute_position: function
## compute_statistic: function
## data: waiver
## draw_geom: function
## finish_statistics: function
## geom: <ggproto object: Class GeomRasterAnn, Geom>
## aesthetics: function
## default_aes: uneval
## draw_group: function
## draw_key: function
## draw_layer: function
## draw_panel: function
## extra_params:
## handle_na: function
## non_missing_aes:
## optional_aes:
## parameters: function
## required_aes:
## setup_data: function
## use_defaults: function
## super: <ggproto object: Class Geom>
## geom_params: list
## inherit.aes: TRUE
## layer_data: function
## map_statistic: function
## mapping: NULL
## position: <ggproto object: Class PositionIdentity, Position>
## compute_layer: function
## compute_panel: function
## required_aes:
## setup_data: function
## setup_params: function
## super: <ggproto object: Class Position>
## print: function
## show.legend: NA
## stat: <ggproto object: Class StatIdentity, Stat>
## aesthetics: function
## compute_group: function
## compute_layer: function
## compute_panel: function
## default_aes: uneval
## extra_params: na.rm
## finish_layer: function
## non_missing_aes:
## parameters: function
## required_aes:
## retransform: TRUE
## setup_data: function
## setup_params: function
## super: <ggproto object: Class Stat>
## stat_params: list
## subset: NULL
## super: <ggproto object: Class Layer>
## ..$ :Classes 'LayerInstance', 'Layer', 'ggproto' <ggproto object: Class LayerInstance, Layer>
## aes_params: list
## compute_aesthetics: function
## compute_geom_1: function
## compute_geom_2: function
## compute_position: function
## compute_statistic: function
## data: data.frame
## draw_geom: function
## finish_statistics: function
## geom: <ggproto object: Class GeomRect, Geom>
## aesthetics: function
## default_aes: uneval
## draw_group: function
## draw_key: function
## draw_layer: function
## draw_panel: function
## extra_params: na.rm
## handle_na: function
## non_missing_aes:
## optional_aes:
## parameters: function
## required_aes: xmin xmax ymin ymax
## setup_data: function
## use_defaults: function
## super: <ggproto object: Class Geom>
## geom_params: list
## inherit.aes: FALSE
## layer_data: function
## map_statistic: function
## mapping: uneval
## position: <ggproto object: Class PositionIdentity, Position>
## compute_layer: function
## compute_panel: function
## required_aes:
## setup_data: function
## setup_params: function
## super: <ggproto object: Class Position>
## print: function
## show.legend: FALSE
## stat: <ggproto object: Class StatIdentity, Stat>
## aesthetics: function
## compute_group: function
## compute_layer: function
## compute_panel: function
## default_aes: uneval
## extra_params: na.rm
## finish_layer: function
## non_missing_aes:
## parameters: function
## required_aes:
## retransform: TRUE
## setup_data: function
## setup_params: function
## super: <ggproto object: Class Stat>
## stat_params: list
## subset: NULL
## super: <ggproto object: Class Layer>
## $ scales :Classes 'ScalesList', 'ggproto' <ggproto object: Class ScalesList>
## add: function
## clone: function
## find: function
## get_scales: function
## has_scale: function
## input: function
## n: function
## non_position_scales: function
## scales: list
## super: <ggproto object: Class ScalesList>
## $ mapping :List of 2
## ..$ x: symbol lon
## ..$ y: symbol lat
## $ theme : list()
## $ coordinates:Classes 'CoordMap', 'Coord', 'ggproto' <ggproto object: Class CoordMap, Coord>
## aspect: function
## distance: function
## is_linear: function
## labels: function
## limits: list
## orientation: NULL
## params: list
## projection: mercator
## range: function
## render_axis_h: function
## render_axis_v: function
## render_bg: function
## render_fg: function
## train: function
## transform: function
## super: <ggproto object: Class CoordMap, Coord>
## $ facet :Classes 'FacetNull', 'Facet', 'ggproto' <ggproto object: Class FacetNull, Facet>
## compute_layout: function
## draw_back: function
## draw_front: function
## draw_labels: function
## draw_panels: function
## finish_data: function
## init_scales: function
## map: function
## map_data: function
## params: list
## render_back: function
## render_front: function
## render_panels: function
## setup_data: function
## setup_params: function
## shrink: TRUE
## train: function
## train_positions: function
## train_scales: function
## vars: function
## super: <ggproto object: Class FacetNull, Facet>
## $ plot_env :<environment: 0x0000000018e0e8b0>
## $ labels :List of 6
## ..$ x : chr "lon"
## ..$ y : chr "lat"
## ..$ xmin: chr "xmin"
## ..$ xmax: chr "xmax"
## ..$ ymin: chr "ymin"
## ..$ ymax: chr "ymax"
## - attr(*, "class")= chr [1:2] "gg" "ggplot"im + geom_path(aes(x=lon, y=lat), data=gpx)